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We perform direct large molecular dynamics simulations of homogeneous SPC/E water nucleation, 
using up to ~ 4 • 10® molecules. Our large system sizes allow us to measure extremely low and 
accurate nucleation rates, down to ~ 10^®cm“®s“^, helping close the gap between experimentally 
measured rates ~ 10^^ cm“®s“'^. We are also able to precisely measure size distributions, sticking 
efficiencies, cluster temperatures, and cluster internal densities. We introduce a new functional form 
to implement the Yasuoka-Matsumoto nucleation rate measurement technique (threshold method). 
Comparison to nucleation models shows that classical nucleation theory over-estimates nucleation 
rates by a few orders of magnitude. The semi-phenomenological nucleation model does better, 
under-predicting rates by at worst, a factor of 24. Unlike what has been observed in Lennard-Jones 
simulations, post-critical clusters have temperatures consistent with the run average temperature. 
Also, we observe that post-critical clusters have densities very slightly higher, ~ 5%, than bulk 
liquid. We re-calibrate a Hale-type J vs. S scaling relation using both experimental and simulation 
data, finding remarkable consistency in over 30 orders of magnitude in the nucleation rate range, 
and 180 K in the temperature range. 

PACS numbers: 05.10.-a, 05.70.Fh, 05.70.Ln, 05.70.Np, 36.40.Ei, 36.40.Qv, 64.60.qe, 64.70.F, 64.60.Kw, 


64.10.+h, 68.35.Md, 83.10.Mj, 83.10.Rs, 83.10.Tv 

I. INTRODUCTION 

The vapor-to-liquid transition of water is a common phe¬ 
nomenon in nature, relevant to many areas of technology 
and science. Attempts to predict the rate of homogeneous 
water nucleation often fail because of the lack of under¬ 
standing of the properties of the tiny seeds of the interme¬ 
diate phase, which are not necessarily large enough to have 
reached the bulk liquid properties. The relevant proper¬ 
ties of the tiny clusters which affect predicted nucleation 
rates include surface tension, temperature, and density. 
Molecular dynamics simulation has proven to be a pow¬ 
erful test of thermodynamic analytical nucleation models, 
now that codes are efficient enough, and computers fast 
enough. Realistic, atmospheric nucleation rates are too 
low to be possible in direct computer simulations, due to 
the large number of molecules required. The lowest water 
nucleation rates performed in simulations and reported in 
the literature are ~ cm“^s’ usually beyond 

the spinodal limit. Laboratory water nucleation rates on 
the other hand are far lower - usually < 10^°cm“^s“^, al¬ 
though a few experiments have managed to measure far 
higher rates ~ 10^^ cm’ Our simulations of ho¬ 

mogeneous SPC/E water nucleation, which we report on in 
this paper, manage to close the gap considerably, resolving 
nucleation rates down to ~ 10^®cm“^s“^. 

Nucleation models, which seek to provide explanations 
and predictions for nucleation rates, have a long history of 
falling short when compared to experimental result s[3l|4l[6l- 
[H]. For the case of water, rate predictions from the classi¬ 
cal nucleation theory disagree with experimental measure¬ 
ments by factors of 10 l-103[5H7l[I^^23]. These models also 
have difficulty when predicting rates measured in numeri¬ 
cal molecular dynamics nucleation simulation experiments. 
However, with molecular simulation, one can make mea¬ 



FIG. 1. A slice through the simulation T325f after 581ns. The 
color-map indicates the density, meaning that the white spots 
represent large clusters. By the end of the simulation, the largest 
cluster in this run has 527 members. This simulation box is 
~ 10pm X 10pm X 10pm, although only a thin slice into the 
2 —direction is visible. 


surements more detailed and accurate than what’s possible 
in laboratory experiments. Size distributions, nucleation 
rates, cluster densities, temperatures, and even cluster pres¬ 
sures, shapes, angular momenta, and surface tension mea¬ 
surements are possible. Understanding the properties of the 
tiny yet complex, many-body clusters which form is vital 
for the development of a complete and successful thermo¬ 
dynamic description of the phase transformation [2¥| . Sim¬ 
ulations allow us to identify the shortcomings in the as- 



sumptions made by existing nucleation models, and sug¬ 
gest ways they may be improved. Cluster properties are 
noisy, necessitating large systems with many millions of 
molecules. This demands costly compute power, and only 
recently have some of these direct measurement techniques 
become possible [2HH27] . 

Direct vapor-to-liquid molecular dynamics simulation for 
a Lennard-Jones fluid has become a popular exercise due 
the computational accessibility of the short-range, single¬ 
site potential[26l l28H^ . Water is significantly more de¬ 
manding. For the same system size, more complicated 
molecular interaction potentials like SPC/E [55H55] and 
TiP4P[aiMiiiz] necessitate a few orders of magnitude more 
computational power than a pure Lennard-Jones simula¬ 
tion. An exception is mW water, a comparatively simple 
monoatomic single-site water mo del [551144) . MD nucleation 
simulations of mW water have been carried out, yet only on 
small systems with relatively high nucleation rates [441145] . 
The monoatomic water model proposed by Zipoli et al. 
(2013)115] offers similar advantages. However, we found 
that short-range potentials require extremely long equili¬ 
bration times to form the correct equilibrium abundance 
of small clusters (dimers, trimers, etc.) in a supersatu¬ 
rated vapor, because interactions are rare, especially the 
three body encounters required for dimer formation. This 
drawback makes it computationally expensive to simulate 
realistic, steady state vapor-to-liquid with such short-range 
potentials - despite their low cost per time-step - and we 
do not use them in this work. Matsubara et al. (2007) [34] 
simulate homogeneous vapor-to-liquid nucleation using the 
SPC/E water model and include a Lennard-Jones carrier 
gas, measuring rates down to 2.3 • 10^®cm“^s“^. SPC/E 
simulations by Tanaka et al. (2014) [1] manage to reach 
nucleation rates 3 • 10^^cm“^s“^. Both efforts addition¬ 
ally measure critical cluster sizes, formation energies, size 
distributions and sticking probabilities for systems in the 
T = 300 — 390 K, providing ample opportunity for model 
comparison and development. 

In this study, we continue in similar spirit, yet simulat¬ 
ing the SPC/E water vapor-to-liquid phase change in even 
larger computational volumes using longer time integra¬ 
tions. This allows for the measurement of lower nucleation 
rates than previously possible by a few orders of magnitude, 
and for the first time, measurements of naturally-formed 
SPC/E cluster density and temperature profiles. Our re¬ 
sults provide opportunities for the verification and calibra¬ 
tion of the standard assumptions which go into nucleation 
models, in a previously unexplored temperature and satu¬ 
ration regime. 


II. SIMULATIONS 

A. Simulation code, setup and parameters 

We use the molecular dynamics SPC/E[53] water model. 
SPC/E is a rigid 3-site model, which registers Coulombic 
interactions, as well as polarization corrections to each site. 


and further adds a Lennard-Jones component to the oxygen 
atom potential. 

The Large-scale Atomic/Molecular Massively Parallel 
Simulator (or LAMMPS) computer program[47|, developed 
at the Sandia National Laboratories and distributed under 
the GPL license, was used to perform the SPC/E simula¬ 
tions. We have verified that our runs produce the same 
results as found in similar, yet smaller SPC/E numeri¬ 
cal nucleation experiments[T|, which used an independent 
molecular dynamics code. We cut the short range Lennard- 
Jones component to the force field off at 9.8 A. Eor these 
forces, as well as the others, the interactions are computed 
directly on per atom. However, the SPC/E Coulombic in¬ 
teractions are long range, and so after 63 A, the spectral 
solver takes over, and the interactions computed in recipro¬ 
cal space. LAMMPS uses a particle-particle/particle-mesh 
solver. The solver maps the atom charges onto a mesh, 
solves the Poisson equation (Maxwell’s equation for the 
electric field) by performing a 3D fast Fourier transform, 
then interpolates the electric fields on the mesh points back 
onto the atom positions |i7M5] . SPC/E molecule rigidity 
is ensured through the use of the SHAKE algorithm|5D]. 
We choose an integration time-step of Af = 2fs, common 
for SPC/E water simulations[T|. A typical simulation runs 
for 72 hours on 1024 cores. Our largest simulation ran for 
1000 hours on 8192 cores on the Piz Daint supercomputer 
at Centro Svizzero di Calculo Scientifico (CSCS), perform¬ 
ing 3 • 10® integration time-steps. 

The simulation box has periodic boundary conditions. 
Initially the molecules are given random non-overlapping 
positions and random velocities. This is done at 1000 K, 
after which the ensemble is cooled and the box size ex¬ 
panded until the simulation reaches the target temperature 
and pressure. The run continues in this state under NVT 
conditions, regulated by a Nose-Hoover thermostat [5TH55] 
with temperature damping timescales of 1000 fs. 

At this stage the gas is allowed to equilibrate for a fixed 
amount of time - dependent on the run temperature (Re¬ 
fer to table |I| for the chosen equilibration timescales at 
each temperature). During this phase the subcritical clus¬ 
ter equilibrium distribution forms. Around this stage we 
begin to make nucleation rate, size distribution, and cluster 
growth rate measurements. Eor most runs, the nucleation 
rate is low enough that unnatural effects from the interven¬ 
tions due to the thermostat are minimal. Our nucleation 
rates are low enough that the latent heat of transformation 
in the simulations is extremely small, resulting in only a 
faint influence from the thermostat. Our largest run sees 
a total energy increase of ^ 0.1% over the steady-state 
phase, i.e. our simulations are very close to NVE (micro- 
canonical) ensembles. 

The first few columns of table|H|lists the runs which were 
carried out, their target temperatures, box sizes, number of 
molecules, and their run times. 


2 


TABLE 1. Thermophysical quantities and parameters at each 
temperature. The vapor equilibrium pressure Pv, the planar 
surface tension 7 , and the bulk liquid density pb for SPC/E 
water are determined from the fitting functions in Matsubara et 
al. (2007) [33]. T] and ^ are nucleation model parameters [3D], te 
the time over which we allow our simulations to equilibrate into 
the steady-state before taking size distribution measurements. 


T 

R 


7 

Pb 

ro 

n 


te 

[K] 

[dyn/ 

cm^] 

[dyn/cm] 

[g/cm®l 

[ 10 “® cm] 



[ns] 

300 

8.9- 

10 ^ 

53.4 

0.997 

1.93 

6.05 

8.95 

25 

325 

4.1 • 

10 "^ 

50.1 

0.982 

1.94 

5.29 

7.47 

20 

350 

1.49- 

10 ® 

46.6 

0.966 

1.95 

4.63 

6.31 

10 

375 

4.47- 

10 ® 

42.9 

0.946 

1.97 

4.02 

5.38 

6 


B. Simulation Analysis 


We use the simple Stillinger criterion]^ (also known as 
the friends-of-friends method) to identify clusters. As the 
simulation runs, the cluster size distribution is regularly 
calculated and outputted, typically resulting in > 1000 
size distribution histograms per simulation. The linking 
length was set at OA for all runs, and was tested to yield 
stable size distributions under convergence tests. Further¬ 
more, this choice yields a monomer-dimer number ratio 
consistent with what is expected from the second virial 
coefficient applied to the SPC/E interaction potential |M]. 
The regularly-outputted size distributions can then be con¬ 
verted into cluster threshold sizes, whose slopes in the 
steady-state regime are the nucleation rates. Refer to sec¬ 
tion jlllj for further details on the nucleation rate analysis, 
as well as the results. From the nucleation rate vs. super¬ 
saturation ratio landscape, we calculate the critical clus¬ 
ter sizes using the first nucleation theorem [551 IShj .The size 
distributions also allow us to follow the growth rate of the 
largest clusters in each simulation, providing a measure¬ 
ment of the monomer-cluster interaction sticking efficiency 
(refer to section VI). 


The measurements of specific cluster properties and how 
they vary with cluster size is crucial to testing assumptions 
used in theoretical nucleation models. However, because 
of the noisy nature of many of these properties, one needs 
many millions of molecules per simulation in order to re¬ 
solve interesting cluster properties. Due to this limitation, 
we perform cluster temperature and cluster density profile 
measurements only for our largest simulation, which con¬ 
tained ^ 4 X 10® molecules. We perform this at the end 
of the simulation, well-within the steady-state nucleation 
regime. This calls for per-atom outputs of velocity and po¬ 
sition information. Sections IVIIII and m detail how the 
density profile and temperature measurements respectively 
are made, and discuss the results. 


III. NUCLEATION RATES 


We use a modified Yasuoka-Matsumoto method [35] 
(threshold method) to measure nucleation rates. In the 
steady-state nucleation regime, the time rate of increase of 
the number of clusters above a certain size N is the nu¬ 
cleation rate. However, the simulations must equilibrate - 
form the sub-critical size distribution - and properly popu¬ 
late it before they reach the steady state nucleation regime. 
How long the simulations take to transition into the steady 
state regime is not known a priori. Nucleation rates esti¬ 
mated from the first nucleation event alone (e.g. mean first 
passage time or survival probability methods) can be orders 
of magnitude smaller than the true steady state nucleation 
rates[S7| [SSj. Thus we use the following method: To the 
size-threshold curves, we fit the following function, which 
is able to capture the transition from the equilibration to 
the steady-state phase. 


N{> i) = J ■ Af {> i,t) + Af {> i,0), (1) 
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where J is the nucleation rate, tp is the lag time, is 
the relaxation timescale and V the volume of the simula¬ 
tion box. This function captures the system’s transition 
from the initial equilibration phase to the intermediate re¬ 
laxation phase as the clusters begin to form, through to 
the steady-state regime. We count clusters above a certain 
post-critical size N, frequently throughout the simulation, 
and fit the count to this curve, allowing J, tp, and tr to 
vary. Visual inspection of this approach is provided in the 
upper left panels of figures]^ for run T375c respectively. 

Our nucleation rate measurements are listed in [IT] Figure 
plots our simulations’ nucleation rates against supersatu¬ 
ration, and includes comparison to earlier results which 
used smaller simulations and were therefore restricted to 
lower nucleation rates. Estimates for the critical cluster 
sizes using the first nucleation theorem, via 
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are included as annotations. Our nucleation rate results 
can be split into two categories: 


• High temperature (325 K, 350 K, 375 K): Runs at 
these temperatures have nucleation rates in the range 
~ 1019-24 cm“®s“^. These runs have generally low 
errors on the nucleation rates. For the higher nucle¬ 
ation rates, there is an error on the supersaturation, 
as the pressure drops significantly due to the large 
number of clusters forming quickly. 

• Low temperature (300 K): Here we measure nucle¬ 
ation rates in the range ~ 10^®“^"* cm“®s“^. These 
runs suffer from extremely long equilibration periods. 
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TABLE II. Run temperature T, supersaturation S as calculated from the run monomer number density, box length L, molecule 
number N, runtime tend, nucleation rate measured from simulation Jmd, critical cluster size i* from the first nucleation theorem, 
Jsp Semi-phenomenological model prediction, Jmcnt modified classical nucleation theory prediction, i*^Q critical size from the AG 
reconstruction. 
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which continue while the initial large, stable clusters 
are already forming. In other words, the sub-critical 
distribution formation timescale tr is longer than the 
nucleation timescale 1/(J • V). This leads to large 
errors in both the nucleation rate measurements and 
the supersaturation measurements. 


IV. RATE COMPARISON WITH ANALYTICAL 
MODELS 


Nucleation models endeavor to describe the phase change 
process in purely thermodynamic terms. The standard ap¬ 
proach tries to find the balance between the Gibbs free 
energy gain and cost due to the creation of volume and 
surface. The classical nucleation theory fCNTl [551 [5MM] 
is the most basic of them all, and forms the basis upon 
which many appendages have since been added. In the 
CNT, the surface energy term in the Gibbs free energy is 
simply calculated using the planar surface tension, with no 
additional corrections. The GNT nucleation rate is [SS] 


•^CNT — 




2567r^rQ7^ 
2nkBTf {log Sf 


(4) 

where m is the molecular mass, 7 the planar surface tension 
at the run temperature, Pg the monomer partial pressure 
in the simulation box (assuming an ideal gas) gas pressure 


and S the supersaturation 


S = 


Ps 

Pv’ 


(5) 


where Pv is the equilibrium vapor pressure at the run tem¬ 
perature. Tq is a characteristic molecular radius: 


ro = 



( 6 ) 


where pi is the bulk liquid density at the run temperature. 
Table |l] includes the thermodynamic variables for SPG/E 
water, which we use in our analysis and comparison to nu¬ 
cleation models. The CNT predictions for the nucleation 
rates of SPC/E water at our runs’ supersaturations are 
shown as solid curves in figure We find that the CNT 
predicts too-high nucleation rates by factors of 10 ^“^. 

Various authors [ill EH employ a 2 -parameter, temper¬ 
ature dependent correction factor 


JcNT exp 



(7) 


The Manka et al. (2010) [H (see their figure 4) lami¬ 
nar flow diffusion chamber experiments find that the pa¬ 
rameter pair {A,B) = (—27.56,6500 K) corresponds to 
a global fit of their results and previous experiments [ 3 - 
O HTHm EH Ell EH]- With these parameters the CNT 
rate prediction remains unchanged at a temperature of 
—BjA = 235.8 K, and they still increase with tempera¬ 
ture (at a fixed S'), but less strongly than in CNT. These 
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FIG. 2. For the 768000 molecule simulation T375c, clockwise from the upper left panel: (1) Nucleation rate measurements, (2) 
largest-cluster growth curve, (3) number density and monomer partial pressure, and (4) concentrations, all over the entire run 
period. Nucleation rates are measured by counting the number of clusters above a specified threshold size, at periodic time intervals. 
The steady-state regime slope is the nucleation rate. (Refer to section IIB and equation 0.) The dotted vertical lines ind icat e the 
lag times for each size. The cluster sticking probabilities a are calculated from the measured slopes, di/dt, using equation | |10[ ). We 
can see that this run took ~ 20 ns to equilibrate, and spent another ~ 5 — 15 ns in the lag phase before reaching the steady-state 
regime (for the chosen threshold sizes of N = 50 and N = 90). The probability that a cluster-monomer encounter results in the 
cluster growing by one molecule is the sticking probability a = 0 . 22 . 


corrected CNT predictions, when extrapolated to our su¬ 
persaturations (dashed curves in figure]^ under-predict our 
measurements by 1-3 orders of magnitude. Using our data 
at temperatures T — 325,350, 375 K to determine the best- 
fit parameter pair, we find {A,B) = (—20.5,6100) K. With 
our parameter pair the CNT rate prediction remains un¬ 
changed at a temperature of —BjA = 297.6 K, and they 
increase with temperature at an rate between CNT and the 
Manka et al. model. However, we note that the resulting 
curves (dotted lines in figure]^ are not quite steep enough - 
casting doubt on whether a purely temperature-dependent 
correction is sufficient in this high supersaturation regime. 

The Modihed Classical Nucleation Theory (MCNT)[29] 
implements a minor modification to the CNT, namely, it 
stipulates that the free energy of formation of a cluster of 
size one is zero. This results in a free energy shift for all 
cluster sizes. Like the CNT, the MCNT over-predicts the 
nucleation rates and here the differences are even slightly 
larger (factor of 10^“"*). Figureshows the ratio between 
the MCNT model predictions (red markers) and the direct 
MD measurements. Refer to table HIl for the MCNT model 
nucleation rate predictions. 

The Semi-Phenomenological model (SP)[SH [5M72] at¬ 
taches a ~ 1/i? (or ^ correction to the surface ten¬ 


sion, where R is the cluster size under the assumption of 
sphericity. This radial dependence is functionally equiva¬ 
lent to that introduced by the Tolman length[TTl [731 173] . 
although the motivation is different: the coefficient to this 
term is set by the second virial coefficient 772^3! that 
the dimer number density is correctly predicted. The nu¬ 
cleation rate predictions for the SP model relative to the 
measured values are plotted with red markers in figure 
The predictions at T = 300 K are somewhat accurate - 
within a factor of 5 of the measured values, although, as 
noted in section in the measurements at these tempera¬ 
tures carry significant uncertainty in the nucleation rate. 
At the higher temperatures, the SP model under-predicts 
the measured MD rates by factors of 4 — 80. Table |TT| lists 
the SP model nucleation rate predictions. 


V. NUCLEATION RATE SCALING 

In this section we examine the scaling of the nucleation 
rates. For the case of water. Hale (2005) [75] (and simi¬ 
larly for Lennard-Jones in Hale (2010) |76|) uses a scaling 
relation[77j for experimentally measured nucleation rates 
over the range J = 10^“^° cm“^s“^ of 
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FIG. 3. The filled solid markers are the nucleation rates we 
measure from our simulations. The solid curves correspond to 
the classical nucleation theory predictions. The dashed curve 
includes the CNT correction factor 0 used in Manka et al. 
(2011)[6], and the dotted one is our best-fit correction factor 
(the fit excludes the runs at T = 300 K). 


FIG. 4. A comparison between the measured molecular dynam¬ 
ics nucleation rates, and those predicted by the analytic SP and 
MCNT nucleation models, ‘new’ and ‘previous’ refer to the 
simulations detailed in this paper, and those by Tanaka et al. 
(2014) [I] respectively. While the SP model is more reliable than 
the MCNT, the SP model shows a trend of rate under-prediction 
for lower supersaturations. 


In 5 

(Tc/T-1)1-5 • 


( 8 ) 


Tanaka et al. (2014) [3T] showed that this scaling relation 
works well for large scale Lennard-Jones simulations and 
Argon laboratory experiments, albeit with an exponent of 
1.3 instead of 1.5. We confirm that the same scaling rela¬ 
tion ([^ applies well to our SPC/E water nucleation rate 
measurements. However, we find that the combined nucle¬ 
ation rates from both SPC/E simulations and laboratory 
experiments with water are even better scaled by 

In'8' 

(Tc/T-1)1-7 • ^ ^ 



Figure shows the nucleation rates as a function of (§. 
This empirical scaling relation seems to work well over a 
surprisingly wide nucleation rate range - from J = 10~^ 
to J = 10^®cm“5g-i for both MD simulations and experi¬ 
ments. The results from the MD simulations join smoothly 
with the experiments with the scaling by In S/{Tc/T— l)i-7. 
Figurej^also shows, using solid curves, the nucleation rates 
predicted by the SP model for various temperatures. This 
scaling relation also works very well for the SP model. 


FIG. 5. Nucleation rates as function of lnS/{T /Tc —l)i’7 for MD 
simulations and experiments [U |31 [S] [HI [HI (131 HU ■ Our simula¬ 
tions are filled circles and the previous simulations ‘-I-’ symbols. 
The solid curves show the SP model for various temperatures 
(210, 315, 350, and 380K). For thermodynamic quantities such 
as the surface tension and the saturated vapor pressure, we use 
those of the SPG/E at 315, 350, and 380 K for the comparison 
with the MD results, while real water at 210 K (pink curve) for 
the comparison with the experiment. 
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VI. STICKING PROBABILITIES 


The sticking probabilities a can be calculated from the 
rate at which large, stable clusters grow. For each run, 
we observe the size of the largest cluster, and measure 
its growth rate di/dt over the second half of the simula¬ 
tion. Early on in the simulations, before stable clusters 
have formed, the largest designation jumps between clus¬ 
ters. However, the first stable cluster to form is likely to re¬ 
main the largest until the end of the simulation. The upper 
right panel of figure shows our cluster growth rate mea¬ 
surements for run T375c. We find the cluster size i (t) to be 
strongly cubic within the steady-state regime. The cluster 
growth rate is therefore proportional to the surface area. 
This is consistent with what has been found in Lennard- 
Jones nucleation simulations [291 [30] . We may determine a 
[T115D] from 


a = 


dTrrguthU (1) 


'-5 


dt 


( 10 ) 


The supersaturation S dependence includes the effect of 
the evaporation of molecules from the clusters into the 
gas. We list the measured sticking probability results in 
table [n] Sticking probability results for our low temper¬ 
ature T = 300 K runs are somewhat unreliable and can 
exceed unity due to the large number of dimers, trimers, 
and tetrames which also contribute to cluster growth. Eq 
(10) considers the accretion and evaporation of monomers 


only. Our sticking probability measurements are consistent 
with those measured at slightly higher supersaturations in 
Tanaka et al. (2014) [T]. The upper panel of figure plots 
a against S. The sticking probability is a necessary pre¬ 
requisite for performing the AG landscape reconstruction 
procedure for post-critical clusters (see section VII). 


While we are, due to computational constraints, unable 
to probe the low nucleation rates observed in laboratory 
experiments, it is possible to measure cluster growth rates 
under laboratory conditions. We have performed an addi¬ 
tional simulation from the end state of T325c, in which we 
measured a sticking probability a = 0.26. We target the 
temperature and saturation conditions found in Brus et al. 
(2008)125]. Using a Nose-Hoover thermostat we maintain 
the temperature, and gently increase the box size until the 
supersaturation S = 2.5, after which we continue running 
for 60 ns. Under these low pressure conditions, no new clus¬ 
ters nucleate (Brus et al. (2008) [22] report nucleation rates 
^ 10^ cm“^s“^) due to our comparatively small and short¬ 
lived system. However, clusters which had previously nu¬ 
cleated and then grown under the original conditions, per¬ 
sist. Using the largest of these still-post critical clusters, 
we measure a decreased growth rate: an slope shal¬ 
lower by a factor of ~ 7. Including this, and the reduced 
(by a factor ~ 3) monomer number density into Eq. (101 
gives a sticking probability for these laboratory-like growth 
rates of a = 0.15. In nucleation models the sticking effi¬ 
ciency is usually taken to be unity, entering linearly in the 
transition growth rate (typically denoted i?^), as a prefac¬ 
tor to the AGi exponent. We find that in the T = 325 K 



1 2 3 

In S 



In S (T/273K)^-^ 


FIG. 6. Upper panel: sticking probability measurements for 
our simulations (solid markers) and from previous simulations 
(‘-I-’ symbols) [1]. Our simulations continue the expected trend 
of lower growth rates with decreasing gas pressure. The lower 
panel shows a temperature-dependent scaling relation which re¬ 
duces the results to a single curve. The blue donut marker 
indicates the sticking probability measured under experimental 
saturations (refer to the end of section VII. 


and S = 2.5 regime, the water monomer-cluster sticking 
efficiency is approximately one seventh of what is usually 
used in model predictions, implying an expected lowering 
of predicted nucleation rates by the same factor. 


VII. FREE ENERGY RECONSTRUCTION 

In this section, we evaluate the formation free energy 
of a cluster AGi{S) directly from our molecular dynamics 
simulations, even for post-critical cluster sizes. We obtain 
AGi (S') from the equilibrium size distribution of the clus¬ 
ter. The equilibrium size distribution can be obtained us¬ 
ing the steady state size distribution, the accretion rate of 
molecule on a cluster, and the nucleation rate, all of which 
can be measured directly in from the MD simulations. Re¬ 
fer to Tanaka et al. (2014) |3T] for a thorough explanation 
of the technique. The cluster size distributions are mea¬ 
sured in the MD simulations and time-averaged over the 
steady state nucleation phase. In the accretion rate, we 
use the value of the sticking probability obtained from MD 
simulations. With the use of them, we reconstruct the full 
equilibrium size distribution (at all sizes i, where we have 
good abundance estimates, including i >> i*) and then the 
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entire free energy function AGi{S = 1). We can further de¬ 
rive AGi{S = 1), which is a surface term corresponding to 
the work required to form the vapor-liquid interface, by 
subtracting the volume term from AGi{S): 

AGi{S = l) = AG,{S) + (i-l)\nS. (11) 

Figure shows AG{S = 1) obtained from the MD results 
at 375 K, and various supersaturations. Since AG{S = 
1 ) is supersaturation independent the values from all runs 
should overlap at all sizes. However, our simulation data is 
only good enough for accurate abundance estimates below 
a certain cluster size, which depends on the run properties. 
The highest nucleation rates run of these (T375a) produced 
a large number of clusters over the entire plotted size range 
and allows the most reliable reconstruction of AGi (S). The 
results from the other runs are only accurate at smaller 
sizes, where they overlap with (T375a). Figurej^also shows 
the surface energy AGi{S = 1) divided by that of the CNT, 
i.e., AGi{S = l)/{rii^/^kT). In the figure, we also show the 
results of the SP model, given by 

AG,(S = l)/{r^e/^kT) = l + (12) 


The simulation results deviate from the SP model at 375 K. 
In Figure]^ we can fit the reconstructed AGi{S = 1) with 

AGi{S = l)/{'qi^/^kT) = 1 -h (13) 


using a fitting parameter A = 0.9 for small clusters. 

Figure shows the ratios between the model and (13) 
(H=0.9, 1.0, 1.0 and 1.5 at 375, 350, 325, and 300 K, re¬ 
spectively) and the MD simulations for two cases: one in 
which a = 1 and the other in which a is set to be value 
obtained directly from simulation. In Figure the predic¬ 
tions from the SP model are also shown for comparison. 
We find the new model agrees with the simulations within 
one order of magnitude for all cases. At 375 K this is no 
surprise, since this data was used to to determine the pa¬ 
rameters of our fitting function for the surface term (13). 
The good agreement at the other temperatures is encourag¬ 
ing and might motivate using (13) also to predict nucleation 
rates at different temperatures and supersaturations. 


VIII. CLUSTER DENSITIES 

It has been shown [75] that for spherical clusters, liquid- 
vapor interface densities are well-approximated by 


p{r) = 


1 


Pc + Pg - [pc - Pg) tanh 2 


r — R 


(14) 


where Pc is the number density within the cluster, pg the gas 
number density, R the interface position, and d its width. 
In each cluster’s center-of-mass frame, we bin the spheri¬ 
cal number density, using a bin size of l.sA. The number 
density profiles for clusters of the same size are used to 


make ensemble averages, to which equation (14) can be 


fit. This method of measuring internal cluster densities 



FIG. 7. Reconstructed Gibbs free energy curves shifted to S' = 1 
(top panel). Runs at T = 375 K at different supersaturations 
were used. To get to S = 1 the GNT volume term was sub¬ 
tracted, the resulting AGi{S = 1) can be interpreted as the 
surface term. The bottom panel shows the AGi{S = 1) divided 
by the surface term from GNT. The surface term from the SP 
model and a simple fitting function (Eq. ( |13| )) are shown with 
dotted and solid lines. 


is robust only for clusters which are large enough to pos¬ 
sess a constant density core. Clusters with i < 20 — 30 
are unlikely to have reached a shape well-describable by 
(14), larger clusters’ density profiles on the other hand are 
well-suited to this functional form. Density profile mea¬ 
surements are noisy, and so particularly large runs with 
many clusters in each size bin are necessary for the ensem¬ 
ble average to provide acceptable accuracy. For this reason 
we perform the density profile measurements in our largest 
simulation, T325f, and we do so at the end of the run. Fig¬ 
ure plots Pc against R for clusters in T325f. We observe 
an over-density for clusters between 4 — 5.5A, however the 
clusters approach the bulk liquid values as they grow, al¬ 
though there seems to be a weak overdensity indication of 
~ 5%. This is in contrast to recent Lennard-Jones nucle¬ 
ation simulations [26] which showed cluster densities signifi¬ 
cantly lower than the bulk liquid values. For post-critically 
sized clusters this was attributable to the increased cluster 
temperatures, due to the residual latent heat which had not 
been efficiently redistributed back into the gas. We surmise 
that the good agreement our internal cluster densities have 
with the bulk liquid values to be due to the fact that they 
are in thermal equilibrium with the surrounding gas. 

In the Lennard-Jones case|26], the lowered densities for 
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FIG. 8. Comparisons of nucleation rate from the MD simula¬ 
tions and several model predictions: the SP model (grey filled 
circles), our new surface term fit (Eq. with a = 1 (filled 

circles) and using the a values measured in the MD simulations 
(open circles). 


clusters with i = i* implied larger surface areas, and there¬ 
fore larger-than-expected surface energies, resulting in an 
increased free energy cost to form a critical cluster, which 
lowered nucleation rates from model predictions. We sus¬ 
pect that nucleation rate predictions are more successful for 
SPC/E water than they are for Lennard-Jones because the 
assumption of small clusters possessing the bulk density for 
i = i* \s more realistic for the case of SPC/E water. 


IX. TEMPERATURES 

We define the temperature of an ensemble of atoms from 
their mean kinetic energy 

2 1 
kT = -(Ekinetic) = ^ X! 

i=l 

Using full per-particle velocity information outputted at the 
end of the simulation, we are able to investigate the cluster 
size dependence of temperature. We find that sub-critical 
clusters are at the run average temperature, as observed 
in similar Lennard-Jones simulations [26) . However, con¬ 
trary to what has been observed in Lennard-Jones nucle¬ 
ation simulations, post-critical SPC/E water clusters pos¬ 
sess temperatures consistent with the run average temper¬ 
ature. Figure plots the ensemble average (at each cluster 
size i) of their temperatures against the density profile in¬ 
terface midpoint R (i.e., the cluster radius). 

The latent heat from condensation has been efficiently 
dissipated back into the gas, leaving the post-critical clus- 
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FIG. 9. Cluster temperatures for T325f. We observe duster 
temperatures consistent with the run target temperature. There 
is no signal of residual latent heat for post-critical clusters, con¬ 
trary to what has been observed in Lennard-Jones nucleation 
simulations |26|. We suspect this may be due to the long-range 
nature of the SPC/E interaction potential, which enables effi¬ 
cient energy exchange between members of the cluster and the 
gas. 



FIG. 10. Cluster densities for T325f, plotted against interface 
positions, both determined from density profile fits to (14l. We 
observe an over-density for clusters between 4 — 5.5A, however 
the clusters approach the bulk values as they grow. We are un¬ 
able to robustly implement the fitting procedure for very small 
clusters < 2.5A as they are not yet large enough to have con¬ 
verged to this shape, nor are they spherical. The orange region 
shows the expected critical cluster size, as estimated from the 
first nucleation theorem. 


ters in thermodynamic equilibrium with their surroundings. 
This finding is consistent with the post-critical clusters den¬ 
sity profile measurements, which finds their densities at the 
expected bulk density. Had the clusters significant latent 
heat retention, their densities would be correspondingly 
lower. We conjecture that the efficient kinetic energy ex¬ 
change is effected by the long-range Coulombic interactions 
- even molecules deep within a cluster may exchange energy 
and angular momentum with members of the gas - resulting 
in kinetic energy equipartition on shorter timescales than 
cluster growth rates. A molecule impinging on a cluster im¬ 
parts heat onto into the cluster-system, yet the heat does 
not linger. This may have implications for non-isothermal 
nucleation models |75j, which include latent heat retention 
in the thermodynamic description of growing droplets. 
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X. CONCLUSIONS 

We have performed molecular dynamics simulations of 
SPC/E water, and significantly closed the nucleation-rate 
gap between simulation and experiment, measuring nucle- 
ation rates low as ~ 10“^® cm“^s“^. This is an hitherto 
unexplored saturation and temperature region for water 
nucleation experiments and water nucleation simulation. 
Nucleation rate results in this new regime will provide mod¬ 
els with further testing comparison opportunities, to com¬ 
plement the already-existing lower nucleation rates from 
experiment, and higher nucleation rates from other simu¬ 
lations. We summarize our most significant contributions 
below. 

• We introduce a new functional form, Eq. Q in or¬ 
der to implement the Yasuoka-Matsumoto nucleation 
rate measurement. This modified version smoothly 
captures the system’s transition between the lag 
phase, relaxation phase, and onto the steady-state 
regime. 

• As expected, the CNT over-estimates nucleation rates 
by a few orders of magnitude. The empirical CNT 
correction factor Q [5S], when using the Manka 
et al.[^ best-fit parameter values (calibrated in the 
low nucleation rate, low saturation regime) under¬ 
estimates our rates by a few orders of magnitude. 
When fitting their proposed correction function to 
our results, we find that the slopes are not steep 
enough. We conclude that this empirical and purely 
temperature-dependent correction factor to the CNT 
is not rich enough to reproduce the qualitative be¬ 
havior we observe in our regime. 

• The MCNT nucleation model continues to over¬ 
predict nucleation rates, by factors of up to 10"*. The 
SP model on the other hand, does somewhat bet¬ 
ter, under-predicting rates at worst by a factor of 
24. Despite these failings, we note that these model 
predictions are significantly more accurate than the 
corresponding predictions for the Lennard-Jones fluid 
vapor-to-liquid nucleation [2HI . 


• Performing a cluster growth rate measurement sim¬ 
ulation under laboratory conditions (those found in 
Brus et al. (2008)[2g) of T = 325K and S = 2.5, we 
measure a sticking probability of a = 0.15. This sug¬ 
gests that in this regime, nucleation rate predictions 
from models should lowered by a factor of seven. 

• We find the cluster size i {t) to be strongly cubic 
within the steady-state regime. The cluster growth 
rate is therefore proportional to the surface area, a 
result new to water nucleation. This is consistent 
with what has been found in Lennard-Jones nucle¬ 
ation simulations [ini 130] ■ 

• Unlike Lennard-Jones nucleation simulations, we find 
that post-critical clusters have temperatures consis¬ 
tent with the simulation average temperature: Grow¬ 
ing clusters are in thermal equilibrium with their sur¬ 
roundings. Latent heat is not retained as the clusters 
grow, it is efficiently dissipated back into the gas. 
We suspect this efficiency is due to the long-range 
Coulombic interactions, not present in the Lennard- 
Jones case. This could have an impact on nucleation 
models which include non-isothermal processes into 
the thermophysical modeling of cluster properties (TO). 

• Post critical clusters have densities consistent with 
what is expected from the bulk liquid. There is 
a possible indication of an over-density for clusters 
around the critical size. This would imply a lower- 
than-expected surface area, which lowers the total 
surface energy, decreasing the free energy cost to form 
a critically-sized cluster and would result in higher- 
than-expected nucleation rates. 

• The scaling relation \nS/{T/Tc — 1)^’^ is remarkably 
successful in reducing the 3-parameter T vs. S vs. J 
surface into a 2-parameter curve. It accurately links 
nucleation rates from simulation and experiment from 
over 30 orders of magnitude in the nucleation rate 
range, and a temperature range of 180 K. 
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